Low temperature phase diagram of condensed ^ara-Hydrogen in two dimensions 
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Extensive Patli Integral Monte Carlo simulations of condensed para-Hydrogen in two dimensions 
at low temperature have been carried out. In the zero temperature limit, the system is a crystal 
at equilibrium, with a triangular lattice structure. No metastable liquid phase is observed, as the 
system remains a solid down to the spinodal density, and breaks down into solid clusters at lower 
densities. The equilibrium crystal is found to melt at a temperature close to 7 K. 



PACS numbers: 67.70+n, 61.30.Hn, 68.08.Bc 67. 

A fluid of para- Hydrogen (P-H2) molecules has long 
been regarded as a potential superfluid, due to the light 
mass and the bosonic character of its constituents^. In 
bulk P-H2, however, superfluidity (SF) is not observed, 
because, unlike Helium, molecular Hydrogen solidifies at 
a temperature (T « 14 K) significantly higher than that 
(~ 4 K) at which phenomena such as Bose Condensation 
and, possibly, SF might occur. This is due the depth of 
the attractive well of the potential between two Hydrogen 
molecules, significantly greater than that between two 
Helium atoms. Several, attempts have been madeSiSi^ 
to supercool bulk liquid P-H2, but no observation of SF 
in the bulk phase of P-H2 has been reported to date. 

Reduction of dimensionality is regarded as a plausi- 
ble avenue to the stabilization of a liquid phase of P-H2 
at temperatures sufficiently low that a superfluid transi- 
tion may be observed. This has been the primary moti- 
vation underlying the experimental investigation of ad- 
sorbed films of P-H2 on different substrates. For ex- 
ample, the phase diagram and structure of monolayer 
P-H2 films adsorbed on graphite have been studied by 
various techniquesiSiLSiS One of the most remarkable as- 
pects is that the melting temperature Tm of a solid p- 
H2 monolayer can be significantly less than bulk P-H21S 
This motivates our interest in the study of the phase di- 
agram of P-H2 in two mathematical dimensions, which 
is still relatively unexplored (a systematic study of the 
one-dimensional phase diagram at zero temperature has 
been recently carried out^^). 

Some information has been provided in previous nu- 
merical work by Gordillo and Ceperley^^, and by Wag- 
ner and CeperleyjiSii^ for example, it is known that, in 
the low temperature limit, the equilibrium phase of the 
system is a triangular solid, with a two dimensional (2D) 
density Oo ~ 0.067 A"^. On the other hand, little is 
known about the possible existence of a metastable liquid 
phase, at densities below 60 (i.e., at negative pressure), 
which one may be able to investigate experimentally by 
"stretching" the equilibrium uniform solid phase. Con- 
ceivably, such a liquid phase ought to turn superfluid at 
sufficiently low temperature. 

Microscopic calculations for condensed P-H2 have re- 
cently focused on realistic models of adsorbed films, both 
on graphitei^ii^, as well as on alkali metal substratesiiSiii 
Here, we present results of a theoretical study of the 
phase diagram of condensed P-H2 in 2D, based on path 
integral Monte Carlo (PIMC) simulations. The tempera- 



ture range explored is between 1 and 8 K; we extrapolate 
the results obtained at low temperature to obtain the 
T = thermodynamic equation of state. 

In agreement with previous calculations, we find that 
the equilibrium phase of the system at T=0 is a tri- 
angular crystal; we estimate the equilibrium density 
6*0=0. 0668±0. 0005 A'^. We also estimate the spinodal 
density 9s, namely the lowest density down to which the 
uniform phase can be stretched, before becoming unsta- 
ble against density fluctuations (at which point it breaks 
down into individual clusters). Our computed value of 
9s is 0.0585±0.0010 A~^. Analysis of our numerical re- 
sults suggests that no metastable liquid phase exists in 
this system, at low T; that is, the system remains a solid 
all the way down to 9s, below which it breaks down into 
solid clusters. We also study the melting of the 2D tri- 
angular solid, and determine its melting temperature at 
approximately 6.8 K. 

Although the PIMC method utilized in this work al- 
lows for the sampling of permutations of particles, which 
is essential in order to reproduce in the simulation any ef- 
fect due to quantum statistics, permutations are not seen 
to occur, in the temperature range explored. This is be- 
cause in the crystal phase, the only one observed here at 
low T, permutations are suppressed (as in most solids) 
by the localization of P-H2 molecules; at higher temper- 
ature, on the other hand, though the crystal melts and 
molecules are less localized, they also behave more clas- 
sically, as their thermal wave length decreases. Consis- 
tently with permutations not being important, i.e. P-H2 
molecules obeying Boltzmann statistics in the tempera- 
ture range explored here, no evidence of SF can be seen»i£ 

Our system of interest is modeled as an ensemble of N 
P-H2 molecules, regarded as point particles and whose 
motion is restricted to two physical dimensions. The 
quantum-mechanical many-body Hamiltonian is the fol- 
lowing: 



N 



2m ^ 



(1) 



i<2 



The system is enclosed in a simulation cell shaped as a 
parallelogram of area A, with periodic boundary condi- 
tions in all directions. The density is 9=N /A. In Eq. 

, m is the |^H2 molecular mass and V is the potential 
describing the interaction between two P-H2 molecules, 
only depending on their relative distance. The Silvera- 



Goldman potentiaP^ was chosen to model the interaction 
V, mostly for consistency with existing, comparable cal- 
culations. However, this potential has also been shown 
to provide an acceptable quantitative description of bulk 
condensed p-H2i^242i 

The PIMC method is a numerical (Quantum Monte 
Carlo) technique that allows one to obtain accurate esti- 
mates of physical averages for quantum many-body sys- 
tems at finite temperature. The only input of a PIMC 
calculation is the many-body Hamiltonian (i.e., the 
potential energy function V) . Because thorough descrip- 
tions of PIMC existfS^ it will not be reviewed here. The 
main technical details of this calculation are illustrated 
in Ref. d 

Most of the results provided in this manuscript per- 
tain to a system of N—GA molecules. We also obtained 
results for a system of 144 particles at the T = equi- 
librium density Oo and near the spinodal density Og. At 
the beginning of the simulation, molecules are arranged 
on a triangular lattice. No significant dependence of the 
estimates on the size of the system can be observed for 
the physical quantities studied here, with the exception 
of the energy (see below) . 
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FIG. 1: Energy per molecule e (in K) vs. coverage 9 (A~^), 
computed by PIMC for a 2D P-H2 film of 64 molecules (dia- 
monds). Estimates shown are obtained at a temperature T—2 
K, but are indistinguishable, within statistical uncertainties, 
from those computed at lower temperatures. Dashed line is a 
polynomial fit to the data. Also shown for comparison (cir- 
cles) are the results of the same calculation carried out by 
Gordillo and Ceperley (Ref. |l3), also using PIMC (at T=l 
K) but on a system of 30 molecules or less. 

Fig. n shows computed values of the energy per p- 
H2 molecule (in K) vs the 2D density (coverage) 9, ex- 
pressed in A~^, for a system of N—64: molecules. Energy 
estimates are found to be nearly independent on tem- 
perature, below T « 3 K. Thus, the results shown in 
Fig. n are essentially ground state estimates. A poly- 
nomial fit of the data yields an equilibrium density 9o 
(corresponding to the minimum of the e{0) curve) of 
0.0668±0.0005 A~^. This result is in agreement with 
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FIG. 2: Typical many-particle configuration generated by 
the PIMC simulation at T=2 K, for a system of 144 P-H2 
molecules at a density 9 = 0.067 A~'^, i.e., near the computed 
equilibrium density. Periodic boundary conditions are used 
in both directions. The arrangement of molecules on a trian- 
gular lattice is clearly seen. Each black "cloud" consists of 
320 points, each representing the position of a P-II2 molecule 
along its path in imaginary time (see, for instance, Ref. |2^ . 



a recent, independent Quantum Monte Carlo calculation 
at zero temperature Also shown in Fig. ^are the es- 
timates obtained by Gordillo and Ceperley, who carried 
out PIMC calculations at low temperature on a system 
of thirty molecules or less. The two calculations agree, 
insofar as locating the equilibrium density; there is a nu- 
merical discrepancy between our calculation and theirs, 
which can be attributed to the difficulty of determining 
quantitatively, on a small size system, the contribution 
to the potential energy associated to the periodic im- 
ages of the system outside the simulation cell. For, the 
H2 intermolecular potential has a long-range attractive 
tail whose overall contribution to the potential energy is 
considerably greater than that obtained, for example, for 
condensed Helium at equilibrium, in spite of the fact that 
the interparticle potentials decay as l/r^ at long distance 
in both cases. All of our energy estimates are obtained 
by computing the above-mentioned contribution to the 
potential energy by setting the value of the pair corre- 
lation function to one outside the simulation cell. On 
comparing results obtained on systems with iV=64 and 
A'^=144 particles, we estimate the systematic error on the 
energy values furnished here, due to the finite size of the 
system, to be less than 0.15 K per particle when A^=64, 
and less than 0.03 K per particle for A^=144. 

Fig. 121 shows a typical many-particle configuration, 
generated by the PIMC simulation, for a system of 144 p- 
H2 molecules at a density 6'=0.067 (i.e., close to ^o), 
at a temperature T—2 K. The arrangement of molecules 
on a triangular lattice is clearly seen. Each fuzzy "cloud" 
represents a P-H2 molecule, the typical size of each cloud 
being a measure of quantum delocalization. There is 
essentially no overlap of clouds associated to different 
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FIG. 3: Chemical potential pu(ff) (in K) versus density Q (A~^) 
computed based on Eq. Q for a 2D p-Hg system, using the 
low temperature energy data shown in Fig. Also shown 
(dotted line) is the polynomial fit of the e(ff) results of Fig. Q 



molecules, which is qualitatively an indication that quan- 
tum exchanges are unimportant, and molecules can be 
regarded as obeying essentially Boltzmann statistics. 

Using the polynomial fit to the energy data shown in 
Fig. ^ we computed the low temperature chemical po- 
tential ii{Q) through 



(2) 



The chemical potential is shown in Fig. |31 (solid line). 
The equilibrium density do is identified by the condition 
/Lt(0o) = e(^o). A second density of interest is the spinodal 
{ds), corresponding to the condition dfi/dO — 0. Based 
on our e{e) data, we obtain 61, = 0.0585 ± 0.0010 A-^. 

The spinodal density is the lowest density down to 
which the uniform phase can be stretched (at negative 
pressure), before becoming unstable against density fluc- 
tuations. At any density lower than Og, the uniform film 
breaks down into "puddles" . We can observe this effect 
directly, by examining many-particle configurations gen- 
erated by a PIMC simulation of a system of 144 P-H2 
molecules; an example is shown in Fig. 01 for 2D P-H2 
at a density 0=0.056 A~^, i.e., slightly lower than 6s- 
Such visual observations strongly suggest that such pud- 
dles are not liquid; rather, the system retains therein its 
triangular crystal structure. 

In their 1997 calculation, Gordillo and Ceperley made 
the suggestion that a metastable liquid phase may exist 
at densities below « 0.059 A^^; however, they did not 
attempt to locate the spinodal density in their study>ii 
As it turns out, 0s as obtained in this work lies precisely 
in correspondence of their proposed location of the solid- 
liquid transition. By direct observation of configurations 
such as those shown in Figs. [21 and ^ as well as by exam- 
ining the pair correlation function and the static struc- 
ture factor, we have confirmed the conclusion of Ref. 0, 




FIG. 4: Same as in Fig. Hbut at a density 9 = 0.056 A"^, 
i.e., below the spinodal density ^3 (see text). 




FIG. 5: Energy e (K) per P-H2 molecule as a function of tem- 
perature (K), ai 6 — 0.067 A~^. Numerical estimates pertain 
to a system of A'^=144 P-H2 molecules. The extrapolated T—0 
value is -23.25±0.05 K. The clear change of slope at T~ 6.8 
K corresponds to the melting of the 2D triangular lattice. 
Dashed line is a guide to the eye. 



namely that the system is a crystal, for 9 > 9s- Because 
the uniform phase breaks down below 9s , we further con- 
clude that no metastable liquid phase of P-H2 exists in 
2D, in the T— >0 limit. The system remains a solid all the 
way down to 9s, and breaks down into solid clusters at 
lower densities. In order to study the melting of the equi- 
librium 2D crystal, we computed the energy per particle 
e(T) as a function of temperature for a system of iV=144 
particles. For simplicity, we have kept the density fixed 
at its T—0 equilibrium value, namely 0.067 A~^. The 
results are shown in Fig. \S[ 

At low T, e(r) follows the expected ~ behavior, 
which is consistent with phonons being the low-lying ex- 
citations of the 2D quantum crystal. The extrapolated 
T—0 energy value is -23.25±0.05 K, which is in satis- 
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factory agreement with the value of -23.4 K recently re- 
ported by Cazorla and Boronat, who carried out a T=0 
calculation. We estimate our uncertainty in the deter- 
mination of the potential energy to amount to less than 
0.03 K per molecule, on a system of this size. There is a 
difference of approximately 0.1 K between the above, ex- 
trapolated energy value, and that obtained on a system 
of 64 particles, the latter being the higher one. 

A sudden change of slope is seen to occur at T,„=6.8 
K, where the specific heat c(T) = de/dT takes on a peak, 
which corresponds to the melting of the 2D crystal. This 
value of the melting temperature is comparable to that 
found in studies of P-H2 films on alkali metal substrateg^^, 
and slightly higher than that of a P-H2 surface (computed 
by PIMC),^'^ and of an adsorbed P-H2 monolayer This 
suggests that melting may occur at a lower temperature 
in three dimensions, as zero-point motion of molecules 
in the direction perpendicular to the substrate promotes 
evaporation. In any case, these melting temperatures are 
still too high to expect that a superfluid transition of p- 
H2 may be observed. 

In conclusion, we have carried out extensive PIMC 
studies of 2D condensed P-II2, determining its low tem- 
perature equation of state, inferring its equilibrium and 
spinodal densities, as well as the melting temperature 



of the equilibrium system. Consistently with previous 
study, we found the system to be a triangular crystal at 
equilibrium; this 2D crystal melts at a temperature of 
approximately 6.8 K. 

We have found no evidence of any metastable liq- 
uid phase at low T. The system is found to remain a 
solid all the way down to the spinodal density, below 
which it breaks down into solid clusters. This result, 
perhaps unexpected, raises some doubts about the long 
term prospects of observing a superfluid phase of P-H2 
upon supercooling the liquid. Observing novel phases of 
P-II2 , including a (superfluid) liquid, may require achiev- 
ing a substantial renormalization of the interaction of 
P-H2 molecules, possibly through their interaction with 
the surface electrons of a metal substrate, or of a nanos- 
tructure. An alternate route might be that suggested in 
Ref. IllL namely stabilizing a liquid at low temperature 
by means of an external potential incommensurate with 
the crystal structure of P-H2 . 
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